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ABSTRACT 

We have investigated evolution of magneto-rotational instability (MRI) in 
protoplanetary disks that have radially non-uniform magnetic field such that 
stable and unstable regions coexist initially, and found that a zone in which 
the disk gas rotates with a super-Keplerian velocity emerges as a result of the 
non-uniformly growing MRI turbulence. We have carried out two-dimensional 
resistive MHD simulations with a shearing box model. We found that if the 
spatially averaged magnetic Reynolds number, which is determined by widths of 
the stable and unstable regions in the initial conditions and values of the resis- 
tivity, is smaller than unity, the original Keplerian shear flow is transformed to 
the quasi-steady flow such that more flattened (rigid-rotation in extreme cases) 
velocity profile emerges locally and the outer part of the profile tends to be super- 
Keplerian. Angular momentum and mass transfer due to temporally generated 
MRI turbulence in the initially unstable region is responsible for the transfor- 
mation. In the local super-Keplerian region, migrations due to aerodynamic gas 
drag and tidal interaction with disk gas are reversed. The simulation setting 
corresponds to the regions near the outer and inner edges of a global MRI dead 
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zone in a disk. Therefore, the outer edge of dead zone, as well as the inner edge, 
would be a favorable site to accumulate dust particles to form planetesimals and 
retain planetary embryos against type I migration. 

Subject headings: accretion disks — instabilities — MHD — planetary systems: 
formation — turbulence 



1. Introduction 



The ubiquity of extrasolar planets strongly suggests that planet formation is a common 
process associated with star formation. However, two serious barriers are recognized in planet 
formation theory: meter-size and type I migration barriers. Terrestrial planets or icy cores for 
gas giants that are embedded in protoplanetary disks tend to lose orbital angular momentum 
and migrate inwa rd through tidal inte raction with the disk gas ("type I migration"). Linear 
calculations (e.g., iTanaka et al.l 120021 ) predict that the planets spiral into the host stars on 
timescales < 10 5 years for the minimum-mass solar nebula model, which is much shorter 
than observationally inferred disk lifetime. This is the type I migration barrier for survival 
of planets with an Earth mass or more. 

The meter-size barrier is the barrier for formation of planetesimals. Since rotation 
speed of disk gas is slightly slower than dust grains or planetesimals due to negative radial 
pressure gradient due to global structure of the gas disk and the motions of meter-sized 
bodies (boulders) are marginally coupled to gas motion through aerodynamical gas drag, 
the meter-sized boulders suffer the fastest inwa rd orbital migration and the times cale to 
spiral into the protostar is ~ a hundred orbits (lAdachilll976l ; IWeidenschiUingi 119771 ). This 
timescale would be too short compared to the sticking timescale for the boulders to form 
planetesimals of 1-10 km radius or more, the motions of which are decoupled from t he gas 
motio n, since the meter-sized boulders are expected to stick together only poorly (IBenz 
2000h . 



One way to bypass the meter-size barrier is formation of clumps through self-gravitational 
instability of a dust layer, which can occur on orbital periods (ISafronovlll969l ; iGoldreich fc Ward 
19731 ). However, local Kelvin-Helmholtz instability due to difference in rotation velocities 
between the dust-rich layer (Keplerian) and an overlaying dust-poor layer (sub-Keplerian) 
would prevent dust grains from settling to the midplane to for m the thin enough layer for 
the self-gravitational instability (jWeidenschilling fc Cuzzilll993l ). Furthermore, if global tur- 
bulence exists in the disk, it also stirs up the dust grains from the midplane. Commonly 
observed strong Ha emission from young stars implies protoplanetary disks are viscously 
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evolving accretion disks with turbulent viscosity. One of the most likely s ources for the tur- 
bulen ce is magnetorotational instability (MRI) in weakly ionized disk gas (IBalbus Hawley 
199lli . 



Although turbulence inhibits formation of the thin dust layer near the midplane, dust 
grains could be trapped into turbulent eddies. It is expected that meter-sized boulders are 
concentrated in anticycloni c eddies fe.g.. lBarge fc Sommeria!ll995l ; IChavanisll2000l ; |Johansen et al 
2004 ; llnaba fc Bargdl2006l ). In particular. iJohansen et al.l (120071 ) carried out a simulation of 
evolution of self-gravitating boulders in an MRI turbulent disk and found that lOOOkm-sized 
clumps are formed in eddies. One of the biggest issues in this model is lifetime of eddies. To 
form clumps, the eddies must persist until boulders are accumulated into the eddies and pro- 
duce highly dust-rich regions. Since rapid formation of highly dust-rich regions is required, 
Johansen et al.l (120071 ) assumed meter-sized boulders that show the most favorable behavior 
in the disks. 

Trapping of dust grains is also possible at a locally high-pressure disk radius, because 
a positive radial pressure gradient induces super-Keplerian gas flow in w hich dust grains 



suffer tail winds, while a normal nega t ive gr adient induces head winds (e.g.. iNakagawa et al. 



19861 ; iKlahr fc Linll2005l ). iRice et al.l (120041 ) demonstrated that meter-sized boulders can be 
accumulated in high-density spiral arms of self-gravitating disks. 

The pressure maximum also exists near the inner boundary of MRI "dead" zone. Inner- 
most disk regions are thermally ionized, so that MRI is active there. In outer regions, X rays 
from host stars and cosmic rays penetrate all the way to the midplane of the disk, and the 
ionization degree may be high enough to activate MRI. On the other hand, in the interme- 
diate regions, only surface layers are MRI active a nd the regions near the midplane may be 
inactive ("dead") (jGammielll996l ; ISano et al.ll2000l ). Assuming steady gas accretion through 
a disk, disk gas surface density jumps up at the inner boundary of the dead zone according 
to the change in effective viscosity due to turbulence, so pressure maximum emerges there. 

The pressure maximum, in other wo rds, the outer boundary of a local sup er-Keplerian 
region also halts inward type I migration (ITanaka et al.ll2002l ; iMasset et al.ll2006l ). If residual 
and/or secondary-generated dust grains maintain the dead zone when accretion of planetary 
embryos proceeds, the inner boundary of the dead zone is a favorable cite to retain planetary 
embryos. However, the inner boundary is generally located well inside 1AU. It cannot provide 
building blocks for terrestrial planets in habitable zones and icy planets nor retain the cores 
to form gas giants around solar-type stars, although it would play an important role in 
architecture of short-period extrasolar planets. 



Kretke fc Lin 



(120071 ). however, pointed out that a local dead zone can appear near the 
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ice line and the ice line can be a spatial barrier for dust migration due to gas drag. Because 
dust grains are the most efficient argents for charge recombination, ionization rates and 
thickness of the active layers rapidly decreases across the ice line due to condensation of icy 
dust grains. Correspondingly, in a range of disk accretion rates, a local dead zone appears 
near the ice line and at the outer boundary of the dead zone large amount of icy grains can 
be accumulated and cores are retained to form gas giants. 



Brauer et al.l (120081 ) performed simulations of dust settling and coagulation with radial 
migration due to gas drag and the ice line effect. They found that the dust to gas ra- 
tio incr eases and formation of planetesimals may be efficient near the ice line. Ilda &: Lin 
(I2008a[) showed that if the ionization rate is order of magnitude larger than that predicted 
by iKretke fc Linl (120071 ) due to dust growth [Kretke Linl (120071 ) assumed that all the dust 
grains have /im sizes], cores stop migration at the ice line and they efficiently form gas giants 
without significant reduction of type I migration speed. Note that in order for the ice line 
effect to actually work, disk accretion rate and dust population must be within some ranges 
of parameters such that the thickness of the active layer is comparable to that of dead zone 
near the ice line, although the ranges are not too restricted. 

We here show through MHD simulations that the pressure maximum associated with 
quasi-steady local super-Keplerian rotation may be created in the MRI marginal regions, 
such as the outer boundary of the global dead zone as well as its inner boundary, without 
requirements of persistent turbulent eddies nor the ice line effect. The accumulation of dust 
grains and retention of planetary embryos at the outer boundary of the dead zone have a 
great importance for formation of terrestrial and jovian planets. MRI is controlled by the 
vertic al component of m agnetic field (B z ) penetrating the disks as well as by ionization degree 
(e.g.. lSano et al.ll2000l ). As shown in later sections, non-uniform temporal MRI turbulence 



that occurs in the marginally stable regions transforms disk gas flow into quasi-steady flow 
that has local rigid-rotation regions. This flow pattern is sustained by non-uniform pressure 
gradient produced by mass transfer associated with the temporal turbulence. In the outer 
regions of the local rigid-rotation regions, gas rotation is super-Keplerian. 

Since dead zo nes may shrink a s dust grains grow and surface areas for charge recom- 
bination decrease (jSano et al.l 120001 ) . such marginal state sweeps from outer disk regions to 
inner regions. Furthermore, the density fluctuations due to MRI turbulence induced by the 
enhanced ionization degree l ead to disruptiv e collisions of small planetesimals and the colli- 
sions reproduce dust grains (jlda et al.ll2008l ). The grains lower the ionization degree so that 
disk gas becomes marginally stable against MRI. Then, grain growth starts again. Thus, 
marginally MRI stable state could be maintained i n significant regio ns at < 10AU until 
disk gas is depleted to some degree (Ilda et al.ll2008l ). lOishi et al.l (120071 ) showed the random 
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torques due to the density fluctuations in the surface active layers affect planetesimals in the 
dead zone near the disk midplane. This effect may result in more continuous dust production 
and help maintenance of the marginally stable state. Thus, it is expected that broad regions 
at < 10AU may once experience such local super-Keplerian motions and accumulate dust 
grains to form planetesimals during disk evolution. 

In section [2j we describe the numerical method and the initial conditions. We use a 
local shearing box with non-uniform B z . In section 3.1, the results with a fiducial model 
are shown. In section 3.2, the dependence of the results on initial settings is presented. 
Summarizing the numerical results, we show in section 3.3 that creation of the rigid rotation 
and super-Keplerian regions is regulated by spatially averaged magnetic Reynolds number. 
Section 4 is conclusion and discussion. 



Numerical Model 



2.1. Shearing box model 



We carry out lo cal two-dimensional MHD simulations of proto planetary disks in the 
shearing box model (jWisdom &: Tremaind Il988l ; lHawley et al.lll995l ). The coordinates are 
centered at ro from a host star and rotating with Keplerian angular velocity at r$ (Qq). The 
radial, azimuthal and vertical directions are x, y and z, respectively. Assuming uniformity 
in the ^/-direction, we simulate 2D flow in the x-z plane. From the flow in the x-z plane, 
the evolution of v y is calculated. We will discuss coordinate sizes in the x-z directions (L x 
and L z ) later. Periodic boundary conditions are applied for the x and z directions. In the 
x direction, Keplerian shear motion in the ^-direction is taken into account in the boundary 
condition. 



2.2. Basic equations 

We adopt compressible resistive magnetohydrodynamic (MHD) equations: 

f)v 1 / R 2 \ 1 

+ (vV)v = — V p + — + — (B.V)B-2fi xv + 3fifcx, (1) 
at p \ 8n J Anp 

f + V-(pv) = 0, (2) 

f)T2 

— = Vx[(vxB)-t?VxB], (3) 
P = c 2 s p, (4) 
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where v, P, p and c s are velocity in the rotating frame, pressure, density and sound speed of 
the gas, respectively. The third and forth terms in r.h.s. of equation ([1]) are the Coriolis force 
and tidal force (the difference between the centrifugal force and gravitational force from the 
central star), in which x is a unit vector in the x-direction. B and rj are magnetic field and 
resistivity. For simplicity, we omit the vertical component of the gravity. We include ohmic 
dissipation in equation ([3]), while we neglect ambipolar diffusion. It s effect is weaker than 
the ohmic dissipation in midplane of inner (< 100AU) disk regions (jJml 119961 ). although in 
the surface layer or disk inner edge, dissipation due to ambipolar diffusion may be important 
dChiang fc Murrav-Clav! 120071 ) . 

We scale length and time by disk scale height (H = c s /VL Q ) and VL$ , respectively. Then, 
the independent parameters in the equations are plasma beta (0) and magnetic Reynolds 
number (P m ) defined respectively by 



2c 2 



2 • 



viz 



(5) 
(6) 



where v& = B/^Airp is Alfven velocity and va z is its ^-component. We assume spatially 
uniform and time-independent c s and rj. 



We solve the resistive M HD equations using the CIP-MOCCT method (jYabe fc Aoki 



199ll ; IStone fc Norman! Il992l ) with grid sizes of 0.01H (for dependence of our results on 
resolution, see Appendix 1). We usually integrate the evolution until t ~ 100£7q . Although 
most of numerical studies on nonlinear stages of MRI have assumed ideal MHD, we consider 
low ionization state with finite resistivity. Several numerical simulations ([Fleming fc Stone 
2000l ; lHawley et al.lll996l ) showed that the finite magnetic resistivity reduces growth rates of 
the MRI. 



2.3. Initial conditions 

We assume that gas density is initially uniform (po), so pressure is also uniform (Pq) 
due to the assumption of isothermal gas (constant c s ). We also assign the initial value of 
plasma beta as (3 — 400 uniformly. According to the large (3, we set initial steady flow as a 
uniform Keplerian shear flow, v y = — (3/2)xf2 . 



The remaining parame ter is R m . This value determines the MRI growth rate (IJin 



19961 ; ISano fc Miyamalll999l ). For R m < 1, short wavelength modes are stabilized by ohmic 
dissipation and the growth rate declines. For fi ^ r~ 3//2 , the the most unstable wavelength 
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is flSano fc Mivamalll999h 



n 2\[2ix . 44 
A m . u . ~ 2tt^- ~ ~ —if. 7 

Since MRI would not occur in the disk if H < A mu ., the stabilization condition is _R m < 

#m,crit ^ 0.44. 

We consider marginally stable state for MRI, in which stable and unstable regions co- 
exist non-uniformly. In this study, we assume constant f], and the non- uniformity of i? m is 
set by that of v& z (equivalently, B z ). Since non-uniform MRI is an essential point for the 
emergence of the super-Keplerian flow, similar results are obtained in the case of non-uniform 
rj. We will show the results with non-uniform rj in a next paper. 

In order to represent the non-uniformity of B z with uniform j3, we set initial Bo such 

that 

B = (0,.Bo sinO, B cos 9) (8) 

with uniform Bq and non-uniform 9. As shown in Figured], we set 9 = (cos^ = 1) in the 
middle region (unstable) with radial size L u and 9 = 85 degrees (cos 9 ~ 0.087) in the side 
regions (stable) with individual radial size L s . The box size in the horizontal direction is given 
by L x = L u + 2L S . Transition zones between 9 = and 85 degrees are set to avoid numerical 
instability. We include non-zero azimuthal magnetic component B y because the plasma beta 
and therefore the magnetic pressure are set to be spatially uniform to establish the initial 
equilibration. The azimuthal component B y is calculated even in these two-dimensional 
simulations on the x-z plane, but the assumption of axisymmetry excludes the occurrence 
of magnetic dynamo on the x — y plane due to reconnection of B y and the results here do 
not change even if B y is set to be zero. In our preliminary three-dimensional simulations, 
we find that the MRI growth due to reconnection of B y in the side areas hardly affects the 
features of azimuthally averaged fluid motion. This is because the essential process of the 
transformation from Keplerian flow to quasi-steady non-uniform rotation flow is temporal 
MRI growth and the stabilization of MRI due to established rigid rotation but not due 
to dynamo of magnetic field (see section 3). The effects of three-dimensional flow will be 
addressed in detail in the next paper. 

The above choice indeed situates both stable and unstable regions in the simulation 
box. Substituting equations (jSJ) and (jHJ) into equation (J6]), we obtain 

^ . dg£ . . , 5 (_^_) co , , m 

With a fiducial value rj = 0.002H 2 Q , MRI is initially activated (R m > -R mjCr it) in the middle 
region with 9 < 9 CTit ~ 65 degrees. The exact dispersion relation (Appendix 2) show that 
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wavelengths of all the growing modes exceed the vertical box size L z for 6 < 6 CI i t ~ 79 
degrees. The vertical size L z is set to be ~ 0.5H that is comparable to the most unstable 
wavelength A m . u . (equation [7]). 

We have carried out 79 runs in total with various rj, L u , and L s . The detailed simulation 
parameters are given in Table HJ In some cases, we adopted larger values of r], in which A m . u . 
is larger (equation [7j). Then, we adopt L z as large as A m . u . with R m ~ 1. In all runs, 
random fluctuations are initially given in the velocity field with a maximum amplitude of 
\5v x \ = 0.001c s . 



3. Results 

3.1. The Fiducial Model 

We first present the detailed results from a fiducial model with 77 = 0.002/f 2 fio) L n = 
1A3H, L s = 4.05H, and L z = 0.5H. Compared with the other models, this model initially 
has the largest stable regions. If the whole region has uniform and strong enough magnetic 
field, MRI turbulence is self-sustaining. However, evolution of turbulence is quite different 
in the case of non-uniform magnetic field that we consider here. MRI turbulence appears 
only tentatively, but velocity field is transformed from Keplerian shear flow to another quasi- 
steady flow. 

Time evolution of the magnetic field on the x-z plane is shown in Figure [2k. In the 
panel of tQo = 21, the field lines start to be stretched only in the initially unstable region. At 
tflo = 27, MRI turbulence develops and the turbulence is transported into the initially stable 
regions. Accordingly, the magnetic perturbations become weaker and eventually vanish 
before they reach the boundary of the computing box (before they fully fill the stable region). 

Figure [2b shows evolution of B z averaged over z. Turbulent diffusion smooths out B z , 
however, its level does not go below the critical value for stabilization of MRI (~ 0.2i?o) in the 
initially unstable region (the middle region) even at tQo = 70. Nevertheless, MRI turbulence 
ceases after tQ > 40. This is because in the middle region, rigid rotation is realized as a 
result of the angular momentum transfer during the MRI turbulence there (Fig. [2b) and the 
free-energy source for MRI (differential rotation) is lost. In the other areas, B z is smaller 
than the critical value and MRI does not grow. The stabilization due to rigid rotation is 
discussed with linear analysis in Appendix 2. 

As shown in Figure [2^, the region of the rigid rotation expands as propagation of 
turbulence. This local rigid rotation is self-sustaining, at least, until the remaining B z is 
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diffused out by ohmic dissipation on timescale ~ L^/i] ~ 10 3 f2 (see discussion in section 4). 
As long as certain level of B z is maintained, deviation from the rigid rotation excites MRI 
turbulence and this transfers angular momentum to bring the velocity profile back to the 
rigid rotation state. 

Figure [2]i is the vertically averaged pressure distribution. The MRI turbulence radially 
transports mass, associated with the angular momentum transfer. Since in the unstable 
region, mass transfer is efficient and gas density is proportional to pressure in our isothermal 
model, pressure is decreased in the unstable region while it is increased in the stable regions 
adjacent to the unstable region. The resultant pressure gradients near the boundaries be- 
tween the stable and unstable regions equilibrates with the Coriolis and tidal forces raised 
by a deviation from Keplerian rotation. In Figure [31 we plot radial components of individual 
terms in the r.h.s. of the equation of motion (equation pQ) near the right boundary between 
the stable and unstable regions (x/H = 0.7, z/H = 0.25). At tfl = 20-40, the amplitude 
of magnetic pressure and tension are increased by the development of MRI turbulence and 
gas pressure gradient is also raised according to mass transfer associated with the turbu- 
lence. After the MRI turbulence ceases, the pressure gradient almost equilibrates with the 
Coriolis and tidal forces. Thus, the rigid rotation with pressure variation caused by the MRI 
turbulence is quasi-steady. 

As stated in section 1, dust grains can be trapped in outer edge of local super-Keplerian 
regions, because dust grains suffer tail wind and migration outward in the super-Keplerian 
regions (see sections 1 and 4) while they suffer head wind in the other sub-Keplerian regions. 
In the result of Figure [21 the super-Keplerian region exists at 0.0 < x/H < 2.0 in the 
quasi-steady state. 

3.2. Dependence on widths of stable and unstable regions 

In the limit of L s — > or L n — > oo, the system would have uniform turbulence and 
would not acquire the local rigid rotation. To derive the condition for establishment of the 
local rigid rotation, we perform runs with different values of L s and L u . 

In the fiducial case, L u = 1AH and L s = 4.0H. We first present the results of a series of 
runs with various values of L s and the fixed L u . In model-sll, model-s055, and model-s005, 
L s = 1.1H. 0.55H, and 0.05H, respectively. The other parameters are the same as those 
in the fiducial model. Figure H] describes the evolution of the magnetic field in model-sll. 
Since the stable region is narrower than the fiducial case, the magnetic perturbations that 
arise in the unstable region propagate to the boundary of the computational box before they 
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are dissipated. The whole region including the initially stable region becomes temporarily 
turbulent. However, the perturbations that reach the boundary are weakened and do not 
have enough momentum to go into the unstable region again. The turbulence ceases after 
tQ > 40 and the quasi-steady rigid rotation is obtained as the result, although the pressure 
contrast is smaller than that in the fiducial model (the minimum pressure is ~ 0.65Poj while 
it is ~ 0.33Po in the fiducial model). 

The evolution of the MRI in model-s055 is shown in Figure The narrower stable 
region allows the magnetic perturbations to come back to the unstable region because of 
the periodic boundary condition. The turbulence lives longer, but it eventually vanishes by 
tQo ~ 70. The quasi-steady velocity profile which is more flattened than the Keplerian is 
established, although it is not as distinctive as the rigid rotation in the fiducial model. 

Figure M shows the result of model-s005. Because of the small stable regions, turbulence 
expanded to the entire region does not cease and uniform turbulence is maintained. Efficient 
angular momentum transfer in the entire region prevents velocity field from having a quasi- 
steady super-Keplerian region. 

So far, we have changed values of L s while L u is constant. We also performed a run 
(model-u34) with L u = 3AH that is enlarged from the fiducial model. The other parameters 
are the same as those in the fiducial model. The results are shown in Figure [7J In the 
enlarged unstable region, the magnetic field is stretched enough for reconnection. After the 
reconnection, long-lived magnetic perturbations are sustained around the boundary between 
unstable and stable regions (the panel at tQo = 58.5). The magnetic perturbations reach 
the computational boundary, because the conversion from magnetic energy to kinetic energy 
during the magnetic reconnection adds horizontal motion to the fluid element. However, 
the magnetic perturbations are not strong enough to pass through the stable region and 
come back to the unstable region. As a result, the evolution is similar to model-sll, except 
that the rigid rotation is established in the regions near the boundary between unstable and 
stable regions. The panels at tQ = 61.5 and tQ = 99.0 show that this flow pattern is also 
quasi- steady. 

These results suggest that dissipation of magnetic perturbations during the passage 
through stable regions plays an essential role in creation of rigid rotation regions. We per- 
formed additional series of runs with a larger value of the resistivity, rj = 0.0028fio-P 2 - The 
result of model- 77O with L u = 1AH and L s = 4.0H (the same as the fiducial model) is shown 
in Figure [8J Because of the faster dissipation, the quasi-steady rigid rotation is established 
more early, although the weaker angular momentum transfer leads to smaller pressure con- 
trast in the quasi-steady state. Even if we reduce L s to ~ 0.5H, the magnetic perturbations 
do not arrive at the boundary. We also performed runs with further larger values of the 
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resistivity, r]/(n H 2 ) = 0.0012,0.0036,0.0044,0.0052 and 0.0060 (Table 1) to confirm that 
the quasi-steady rigid rotation is always established for these values of r/. We discuss the 
dependence on the resistivity in the next subsection. 



3.3. Condition for the rigid rotation 



As we have shown, establishment of the quasi-steady rigid rotation depends on L u , L B , 
and rj. Through the results with various L u , L s , and rj as listed in Table 1, we found that 
the results are sum marized by a single parameter, a spatially ave raged magnetic Reynolds 
number, defined by (jSano fc Miyamalll999l ; ISano &: Inutsukal 120041 ) 



Az,ave 



where v A 2 ,ave is evaluated by a spatially averaged vertical magnetic field (B z 
stage. In our simulation setting, B z ^ vc ~ B (L n + 2cos85°L s )/(L u + 2L S ). 



(10) 



at the initial 



In many cases, the initial distributions of vertical magnetic field are spatially smoothed 
out by turbulent diffusion after ~ 10 orbits. If the remaining magnetic field (~ B ZtStve ) 
is still large enough to globally cause MRI turbulence, the transition to the quasi-steady 
state does not occur, as seen in the results of model-s005. Linear theory shows that the 
critical magnetic Reyn olds number for occurrence of MRI is R m ~ 1 (jSano fc Miyamalll999l ; 
Sano fc Inutsukal 120041 ) . Thus, -R miav e regulates the establishment of the quasi-steady rigid 
rotation in our simulations. 

We summarize the results of 91 runs with different parameters in Figure [9] and found 
that evolution of the magnetic field is indeed classified into four types by the values of Rm,a,ve 
as follows: 



Type A (-R mj ave < 0.1): Local MRI turbulence generated in the initially unstable region 
propagates both inward and outward, but the magnetic perturbations are dissipated 
before they reach the boundary of the simulation box. After a few tens of orbits, the 
turbulence vanishes in the entire region and the quasi-steady flow is established. In 
the initially unstable region, rigid rotation flow is resulted in by angular momentum 
transfer due to the MRI turbulence. This class is represented by filled circles in Figure[9] 
and it includes the fiducial model and model-z^O. 

Type B (0.1 < -R m ,a ve ^ 0.5): The magnetic perturbations reach the boundary, but they 
do not intrude back into the original unstable region. The quasi-steady rigid-rotation 
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region appears as in Type A. This class is represented by triangles in Figure [9] and it 
includes model-sll and model-u34. 

Type C (0.5 < R m ^ ve < 1.0): The magnetic perturbations intrude the unstable region after 
the passage through the stable regions. However, the diffused B z is not large enough to 
globally maintain turbulence and the quasi-steady rigid rotation region is still formed, 
although their locations are not necessarily the same as in Type A and B. This class 
is represented by daggers in Figure [9] and it includes model-s055. 

Type D (1.0 < i? mi ave) : Even after the turbulent diffusion, B z is able to maintain the tur- 
bulence in the entire region. Because of the uniform turbulent state, the quasi-steady 
rigid rotation state is not generated. This class is represented by crosses in Figure [9] 
and it includes model-s005. 



4. Conclusion and Discussion 

We have investigated evolution of patchy magneto-rotational instability (MRI) due to 
radially non-uniform magnetic field and found that, under some conditions, the original 
Keplerian shear flow is transformed into quasi-steady profile involving a local rigid-rotation 
regions. The outer parts of the rigid-rotation regions are generally super-Keplerian. Such a 
situation would arise in the outer boundary of MRI dead zone as well as the inner boundary 
in a protoplanetary disk, as discussed in section 1. 

Assuming uniformity in the azimuthal direction of disks, we have carried out two- 
dimensional resistive MHD simulations in a shearing box model with periodic boundary 
conditions. We set up both stable and unstable regions in the box, changing direction of 
the vertical seed magnetic field (B z ) non- uniformly. In the initially unstable region, MRI 
turbulence is generated locally and magnetic perturbations propagate both radially inward 
and outward by the turbulent diffusion. If the unstable region is sufficiently large compared 
with the stable region, the turbulence eventually covers the entire region and the initial 
non-uniformity vanishes. However, if the stable region is relatively large, diffused magnetic 
perturbations no more maintain MRI turbulence. After the turbulence ceases, the initial 
flow of uniform Keplerian shear is transformed into a different quasi-steady state. In the 
quasi-steady state, rigid-rotation is established locally. The deviation from Keplerian shear 
motion is supported by pressure gradient that has been produced also by mass transport 
associated with the tentative turbulence. Through simulations with various initial condi- 
tions, we found that the quasi-steady rigid rotation is established if the spatially averaged 
magnetic Reynolds number satisfies -R m ,ave < 1 in the initial state. 
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Because the center of the local rigid rotation is often Keplerian, super-Keplerian flow 
appears in the outer parts of the rigid rotation region. As explained in section 1, dust grains 
and planetary embryos can be trapped in the boundary between regions of sub- and super- 
Keplerian motion through radial migration induced by aerodynamic gas drag and type I 
migration. The boundary is coincident with pressure maxima in the quasi-steady state. 

The effect of global pressure gradient is included by shifting the initial gas velocity from 
the pure Keplerian speed to slight sub-Kepler. The sift, which is the velocity diff erence 



c\ d\nP 9 / r \3/4 



between the disk gas and the dust grains, is (e.g., lAdachilll976l ; IWeidenschillingi 119770 

- 2 dh\P / /■ \ :! 1 

K 

where d In P/d In r is the global pressure gradient and the temperature distribution in the 
limit of optically thin disks around solar-luminosity stars, T = 280(r/lAU)~ 1//2 , is assumed. 
Since maximum values of Av y in the super-Keplerian regions are > 0Ac s in our results, 
the MRI effects can easily surpass the global pressure gradient effects and super-Keplerian 
regions can emerge even when the initial down-shift in the gas rotation velocity is present. 

If the super-Keplerian region is sustained long enough for dus t grains to accumu- 



late, planetesimals ca n be formed through self-gravitational instability (jYoudin &: Shul 12002 



Johansen et al.ll2006l ). In the case of R m ^ ve < 1 in which MRI turbulence ceases in the 
entire region after a few tens of orbits, we found that B z is still large enough in the rigid- 
rotation region. In that region, MRI is suppressed by disappearance of shear motion but not 
by dissipation of B z (in the regions other than the rigid-rotation region, diffused-out B z is 
smaller than the value for MRI to occur). If the rigid rotation tries to go back to the original 
Keplerian shear motion, MRI turbulence again occurs and it transfers angular momentum to 
recover the rigid rotation. Thus, the rigid rotation and hence the associated super-Keplerian 
rotation are self-sustaining. When the remaining B z is diffused out by ohmic dissipation on 
timescale ~ L^/v ~ lO 3 ^ 1 , such stabilization mechanism is no more effective. Then, the 
rigid rotation can go back to the original Keplerian shear motion by the residual uniform 
viscosity. However, since MRI no more occurs, the residual viscosity would be very small. 
Thus, it is expected that the super-Keplerian regions would survive long enough for accumu- 
lation of dust grains and formation of planetesimals. We also did a calculation starting from 
the end result of the super-Keplerian rotation state, artificially modifying B z to uniform 
distribution. However, we do not see any relaxation of the velocity field back to Keplerian 
rotation within the timescales of 40 orbits. 

In the next paper, we will demonstrate the accumulation of dust grains and discuss the 
effect of the azimuthal magnetic field especially in stable region by three-dimensional simu- 
lation. We will also show the results of non-uniform resistivity case with constant B z , which 
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may be more likely to occur at the outer edge of a dead zone. We find similar appearance 
of super-Keplerian regions, because intrinsic physics to transform initial Keplerian flow to 
quasi-steady non-uniform rotation flow is temporal generation of MRI and the stabilization 
of the MRI due to the established rigid rotation but not due to dynamo of magnetic field 
(see section 3). 

The appearanc e of the super-Keplerian region also halts inward type I migration of 



planetary embryos (jTanaka et al.ll2002l ; iMasset et al.ll2006l ). Since the process we found also 
works at the outer boundary of a dead zone and the outer boundary migrate from ~ 10AU 
to the proximity of the central star, this process may also help the formation of cores massive 
enough to onset runaway gas accretion and retain terrestrial planets against type I migration. 
This may play an important role in frequency of extrasolar gas gia nts and habitable plane ts. 
We will address this issue with sequential planet formation model Jlda fc Linlbopj . l2008al lb}) 
in a future paper. 



We thank detailed helpful comments by an anonymous referee. 



Appendix 1. Resolution test 



We have investigated the effects of numerical resolution on our results using our fiducial 
model. Four cases are studied. The distribution of angular velocity and the average Maxwell 
stress are shown in Figure [10J The most important issue in our results is the emergence of 
the quasi-steady state in which sub- and super-Keplerian areas exist. Figure [TUk shows this 
to be seen even in the worst resolution case. In addition, the Maxwell stress in Figure [TUb 
appears to be converging at resolutions above dx = 0.01H while this quantity tends to 
increase with resolution. 



Fromang fc Papaloizoul (120071 ) showed the decrease of Maxwell stress with increasing 
resolution and stated that the MRI turbulence activity is notoriously ill-behaved in high 
resolution calculation. Our resolution test, however, shows the convergence with increasing 
resolution. This could be because the ohmic dissip ation (resistivity) is kept const ant in our 
resistive MHD simulations while the simulation of iFromang fc Papaloizoul (120071 ) included 
only the numerical resistivity and it deacreases with increasing resolution. However, more 
detailed study is needed to clarify the dif f erence in convergence between our simulation and 
the simulation of IFromang fc Papaloizoul (120071 ) . which is left to our future study. 



This fact will be important when we investigate the motion of dust particles. The degree 
of particle concentration may be depend on the turbulence activity. While the details are 
left for a future study, one may reasonably expect the possibility of dust concentration at 
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the outer-edge of super-Keplerian area. 

Figure [TUb also shows the eligibility of the integration time. The saturation level varies 
only slightly in tfl > 60.0 in the two high resolution cases. The quasi-steady state has 
already been created by this time. These indicate that our choice the integration time 
m < 100.0 is validated. 



Appendix 2. Dispersion relations 



Our simulations show that the MRI stabilization in strong magnetic field with nearly 
rigid rotation, whi ch is consistent wi t h the linear analyses. The linear analysis using ideal 
MHD equations in iBalbus fc Hawleyi (I199ll ) gave the critical wavelength (equation (2.14b) 
in the paper, with neglect of the Brunt-Vaisala frequency): 



|A 2 



crit 



VAz 



dtt 2 



cHnr 



-1/2 



VAz 



(2|g|) 



-1/2 



where q is defined as Q(r) oc r~ q . The perturbations with wavelength shorter than A 2jCr it are 
stable. When the rotation becomes rigid rotation (q — > 0), A 2jCr it becomes large. If A 2jCr it 
is larger than the scale height of a disk, the system is stable, irrespective of magnetic field 
strength. 



With the effect of ohm ic dissipation, the dispersion relation was obtained by iJinl (119961 ) 
and ISano fc Miyamal (Il999l ) as 



where a is a growth rate in units of orbital frequency, q 
an epicyclic frequency defined by 

2 2VL d (r 2 VL) 



" 4 + ^ 



0, (2) 



k z vAz/Q, £ = k 2 r]/Q, and k is 



K 



dr 



(2 - q) 2VL 1 . 



(3) 



This dispersion relation is derived with the assumption of uniform density and negligible 
Brunt-Vaisala frequency. Since the density is almost uniform in the rigid rotation in our 
simulation, we apply this dispersion relation to calculate the predicted a with the quantities 
obtained by our simulation. We plot the temporally and vertically averaged q (= —dv x /dx), 
Maxwell stress (—(B x B y )/A'n-Po) and the evaluated growth rate a in Figure [TTJ The growth 
rate is very small in the middle region of nearly rigid rotation (q 1), although Maxwell 
stress is not small there. Thus, we conclude that MRI is suppressed by established nearly 
rigid rotation, but not by dissipation of magnetic field. 
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Fig. 1. — The configuration of the initial magnetic field. 9 represents the angle from the 
vertical direction (z) to the azimuthal one (y). It radially varies as shown in the right 
panel. With constant B , the initial magnetic field is given by B = (0, B sin 9, B cos9). 
For 7] = 0.002/J 2 r2 , linear calculations suggest that MRI occurs for 9 < 9 cvlt ~ 79 degrees. 
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Fig. 2. — Results of the fiducial model, (a) Time evolution of the magnetic field (black lines) 
and angular velocity v y (contours) on the x-z plane. The maximum (brightest) or minimum 
(darkest) of the tone on the contours are set at the unstable and stable region boundaries. 
Time evolution of vertically averaged values of (b) the vertical magnetic component B z , (c) 
angular velocity v y , and (d) pressure P, as functions of x. In panels (b), (c) and (d), thin 
solid, dotted, dashed and bold lines express the snapshots at tfto = 0.0,27.0,40.0 and 70.0, 
respectively. The unstable region is initially set between the two vertical dotted-lines. 
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Fig. 3. — Time variation of forces exerted on fluid at x/H = 0.7 and z/H = 0.25, at which 
super-Keplerian flow is established. The dotted, dash-dotted, thin dashed and dashed lines 
express gas radial components of pressure gradient, magnetic pressure gradient, magnetic 
tension, sum of gravity and Coriolis forces in the r.h.s. of the equation of motion (equa- 
tion pQ), respectively. The thick solid line is total force. The unit of the forces is HQ^. The 
pressure gradient and the gravity /Coriolis forces are dominated and approximately equili- 
brate with each other. 
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Fig. 4. — Results of model-sll. (a) Time evolution of the magnetic field (solid lines) and 
angular velocity v y (contours) on the x-z plane, (b) Time evolution of vertically averaged 
angular velocity (v y ). The bold, dashed and thin solid lines express the snapshots at tQo — 
0.0,40.0 and 70.0, respectively. The meanings of lines are the same as Fig. [2j 
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Fig. 5. — The same plots as Fig. [4] except for model-s055. 
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Fig. 6. — The same plots as Fig. Hh except for model-s005. 



-23- 




Fig. 8. — The same plots as Fig. H] except for model-r/O. 
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Fig. 9. — The classification of results by the magnetic Reynolds number. Runs with (a) 77 = 
0.0020fi # 2 and (b)r) = 0.0028fi # 2 on the L u -L s plane and (c)all the runs on the r}-L s /L u 
plane. The four types of results, A, B, C and D, are represented by filled circles, triangles, 
daggers and crosses. The solid, dashed and dotted lines express R m ^ vc = 0.1,0.5 and 1.0 
respectively. The runs with numbers express (l)fiducial model; (2)model-sll; (3)model-s054; 
f4 N )model-s005 f5)model-u34; (filmodel-nO. 



-25 - 




Fig. 10. — Results for various resolutions. (a)Close-up snapshots of vertically averaged 
angular velocity^) at t£l = 150.0. (b)Time evolution of the volume-averaged Maxwell 
stress which is normalized by the initial pressure P . Resolutions corresponding to dx = 
dz = 0.005#, 0.01H, 0M5QH and 0.025// are represented by dashed, bold, dash-dotted 
and dotted lines respectively. 
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Fig. 11. — Temporally and vertically averaged shear rate q = —dv x /dt {top), Maxwell stress 
normalized by initial pressure — (B x B y ) /AtiPq (middle) and the growth rate estimated by the 
dispersion relation (equation (2) in Appendix 2) with the simulated values (bottom) of the 
fiducial model. Bold lines are values temporally averaged over t£l = 45.0 — 70.0 and dashed 
lines mean the initial values. 



-27- 



r]/Q H 2 


T / TT 

LJH 


T 1 TT 

LJH 


-^m,ave 


1 i 

result 


AVy/C S 


0.0012 


1.5 


6.5, 6.0, 5.0, 4.0 


0.10, 0.11, 0.13, 0.16 


A, B, B, B 


0.82, 1.2, 0.91, 0.86 






1.0, 0.51 


0.61, 1.1 


C, C 


0.66, 0.51 




2.0 


7.5, 6.0, 4.0 


0.13, 0.16, 0.23 


B, B, B 


1.7, 1.1, 0.98 






1.0, 0.51 


0.89, 1.4 


C, D 


0.68, ... 




2.5 


0.51 


1.7 


D 




0.002 


1 r\ 

1.0 


11 f\ V A 

1.1, 0.54 


n no r\ ,1 o 

0.23, 0.42 


B, B 


0.40, 0.38 




1.2 


3.1 


0.10 


A 


0.57 




1.4 


4.0 1 , 3.1, 1.1 , 0.55 d 


0.096, 0.13, 0.37, 0.64 


A, B, B, C 


0.73, 0.55, 0.48, 0.41 






0.15, 0.10, 0.05 4 


1.2, 1.4, 1.5 


D, D, D 






1.9 


6.0, 4.0, 3.0 


r\ c\r\ A f~\ -1 A f~~\ -1 /"v 

0.094, 0.14, 0.19 


A, B, B 


0.97,0.88, 0.54 






2.0, 1.0, 0.35 


0.28, 0.53, 1.1 


B, C, C 


0.87, 0.72, 0.43 




2.4 


2.0, 0.35 


0.38, 1.3 


B, D 


1.2, ... 




2.9 


10.0, 8.0, 1.0 


0.092, 0.11, 0.81 


A, B, C 


1.2, 1.4, 0.99 




3.4 


12.0, 8.0, 4.0 5 , 0.24 


0.090, 0.14, 0.28, 1.8 


A, B, B, D 


1.0, 1.2, 1.1, ... 




3.9 


3.1 


0.44 


B 


0.77 


0.0028 


1.4 


5.1, 4.0 6 , 2.1, 0.57 


0.15, 0.069, 0.14, 0.45 


A, A, B, C 


0.35, 0.45, 0.28, 0.23 






0.37, 0.17, 0.11 


0.61, 0.87, 1.00 


C, C, D 


0.34, 0.14, ... 




1.9 


A /~\ i \ -i ' ~\ -i A A 

4.0, 3.1, 2.1, 1.1 


0.099, 0.14, 0.20, 0.38 


A, B, B, B 


0.48, 0.31, 0.49, 0.42 






0.57, 0.17, 0.11 


0.61, 1.05, 1.2 


C, C, D 


0.40, 0.25, ... 




2.4 


6.1, 5.1, 4.0, 1.1 


0.088, 0.11, 0.13, 0.49 


A, B, B, B 


0.50, 0.62, 0.66, 0.53 






0.57, 0.37, 0.11 


0.75, 0.92, 1.25 


C, C, D 


0.46, 0.28, ... 




o n 

z.y 


f.l, 6.1, 0.1, 2.1 


n nno nil nio n yi7 
0.093, 0.11, 0.12, 0.4r 


A "D T> "D 
A, D, D, D 


0.60, 0.4Z, 0.50, 0.4o 






1.1, 0.57, 0.27, 0.17 


0.58, 0.85, 1.13, 1.26 


C, C, D, D 


0.56, 0.37, ... 




3.4 


8.1, 7.1, 6.1, 2.1 


0.0097, 0.11, 0.13, 0.39 


A, B, B, B 


1.00, 0.79, 0.58, 0.78 






1.1, 0.37, 0.27 


0.66, 1.1, 1.3 


C, C, D 


0.48, 0.61, ... 




3.9 


8.1, 2.1, 1.1 


0.11, 0.45, 0.73 


B, B, B 


0.46, 0.43, 0.49 


0.0036 


1.8 


4.1 


0.077 


A 


0.28 


0.0044 


1.7 


4.1, 3.1 


0.063, 0.085 


A, A 


0.17, 0.18 






2.1, 1.1, 0.43 


0.13, 0.24, 0.50 


B, B, B 


0.20, 0.18, 0.18 


0.0052 


2.3 


1.3 


0.24 


B 


0.15 


0.0060 


2.4 


3.3, 2.3, 1.3, 0.32 


0.084, 0.36, 0.21, 0.50 


A, B, B, B 


0.22, 0.17, 0.21, 0.20 



-28- 



Table 1: Simulation parameters and results for 91 runs. r], L u , L s , i? mjave , and Av are 
magnetic resistivity, radial width of unstable region, that of stable regions, initial averaged 
magnetic Reynolds number and the maximum deviation from Kepler velocity, respectively. 
For highly turbulent cases (marked "D"), Av y is omitted. The fifth column "result" indicates 
classification of the results (§3.3). Multiple values in the columns correspond to different 
runs. (1) fiducial model; (2) model-sll; (3) model-s055; (4) model-s005; (5) model-u34; (6) 
model-r/O. 



